# This file creates the variables for analyses with ties removed
# Set working directory
setwd("/Users/ericguntermann/Documents/Papers/Representation of Party Preferences/Replication")

set.seed(123)

load("cses.Rda")
cses <- cses1_4
cses <- cses[-which(cses$countryyear %in% c("Japan 2013")),] # Remove upper house election
## Only keep parliamentary and semi-presidential regimes
cses <- cses[!is.na(cses$regime), ]
cses <- cses[cses$regime<2,]
cses <- cses[order(cses$countryyear),]



## Create a matrix of like dislike scores with weights
cses.ld <- cses[,c(4,11:19, 62:64)]

# Remove people who expressed no preferences
cses.ld[is.na(cses.ld)] <- 11
rowtot <- rowSums(cses.ld[,2:10], na.rm=T)
cses.ld <- cses.ld[!rowtot==99,]
cses.ld[cses.ld==11] <-NA
cses <- cses[!rowtot==99,]

## Create overall weight (sample, demographic, and political)
cses.ld$weight <- cses.ld$sampleweight*cses.ld$demoweight*cses.ld$polweight

# Criterion 1: What proportion of people in each country's first choice party ended up in government?
cses.ld2 <- data.frame(cses.ld)
cses.ld2[is.na(cses.ld2)] <- 0

# Create a vector of first choice parties for all individuals
# Remove all respondents with tied first choices

torm <- vector()
for(i in 1:nrow(cses.ld2)){
  maxtie <- max(cses.ld2[i,2:10])
  n <- sum(cses.ld2[i,]==maxtie)
  if(n>1){
    torm[i] <- TRUE
  }else{
    torm[i] <- FALSE
  }
}
mean(torm)
cses.ld2 <- cses.ld2[torm==FALSE,]
# Create a vector of first choice parties for all individuals
library(nnet)
first <- apply(cses.ld2[,2:10],1, which.is.max)

# Cabinet data
incabinet <- cses[,c(4,21:29)]

# Remove rows for respondents with ties

incabinet <- incabinet[torm==FALSE,]

# Is the most liked party in government (for each respondent)?
firstincabinet <- vector()
for(i in 1:length(first)){
firstincabinet[i] <- incabinet[i,first[i]+1]
}

# Get proportion whose most liked party is in cabinet for each election
propfic <- vector()
for (i in 1:length(unique(cses.ld2$countryyear))){
  propfic[i] <- weighted.mean(as.numeric(firstincabinet[cses.ld2$countryyear == unique(cses.ld2$countryyear)[i]]), cses.ld2$weight[cses.ld2$countryyear==unique(cses.ld2$countryyear)[i]], na.rm=T)
}

# Calculate number of parties in government
## Create aggregate cabinet data matrix

incabinetag <- data.frame(matrix(nrow=length(unique(cses$countryyear)), ncol=9))
for (i in 1:length(unique(incabinet$countryyear))){
  for (j in 1:9){
    incabinetag[i,j] <- incabinet[c(incabinet$countryyear==unique(incabinet$countryyear)[i]),j+1][1]
  }
}
colnames(incabinetag) <- c("incabinet1", "incabinet2", "incabinet3", "incabinet4", "incabinet5", "incabinet6", "incabinet7", "incabinet8", "incabinet9")


ningov <- function(x){
  return(sum(x==1,na.rm=T))
}
ngov <- apply(incabinetag[,1:9], 1, ningov)


# Create dataset
criteria <- data.frame(countryyear=unique(cses$countryyear), propfic)
for (i in 1:length(unique(cses$countryyear))){
  criteria$proportional[i] <- cses$proportional[cses$countryyear==unique(cses$countryyear)[i]][1]
}


for (i in 1:length(unique(cses$countryyear))){
  criteria$mdm[i] <-  cses$mdm[cses$countryyear==unique(cses$countryyear)[i]][1]
}
for (i in 1:length(unique(cses$countryyear))){
  criteria$gallagher[i] <- cses$gallagher[cses$countryyear==unique(cses$countryyear)[i]][1]
}

for (i in 1:length(unique(cses$countryyear))){
  criteria$freedomhouse[i] <-  cses$freedomhouse[cses$countryyear==unique(cses$countryyear)[i]][1]
}

for (i in 1:length(unique(cses$countryyear))){
  criteria$gdppercap[i] <-  cses$gdppercap[cses$countryyear==unique(cses$countryyear)[i]][1]
}


criteria$ngov <- ngov
criteria$coalition <- as.factor(ngov>1)
# Adjust number of parties in Australia. The Liberal and National Parties ran together.
criteria$ngov[criteria$countryyear=="Australia 1996"] <- 1
criteria$ngov[criteria$countryyear=="Australia 2004"] <- 1
write.dta(criteria, "criteria_noties.dta")

# Table 8 of Online Appendix

dat <- read.dta("criteria_noties.dta")
criteria <- dat
criteria2 <- dat[-which(dat$countryyear %in% c("Thailand 2001", "Thailand 2011","Japan 1996", "Japan 2013")),]

#Criterion 1

tab1 <- data.frame(round(mean(criteria2$propfic[criteria2$proportional==0])*100,1), round(mean(criteria2$propfic[criteria2$proportional==1])*100,1), round(mean(criteria$propfic)*100,1))
colnames(tab1) <- c("Non-PR", "PR", "Overall")
rownames(tab1) <- c("Proportion most liked in cabinet (%)")


library(xtable)
xtab <- xtable(tab1, caption="Performance of PR and non-PR Systems Criterion 1 (ties removed)")
sink("tabnoties.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), only.contents=T)
sink(NULL)

# Now run models with no ties (tables 9 and 10)

set.seed(123)
library(foreign)

## Criterion 1 (proportional)
criteria <- read.dta("criteria_noties.dta")
criteria <- criteria[-which(criteria$countryyear %in% c("Thailand 2001", "Thailand 2011", "Japan 1996", "Japan 2013")),]
d1 <- with(criteria, list(propfic=propfic, proportional=proportional, gdppercap=gdppercap, freedomhouse=freedomhouse))
d1$N <- length(unique(criteria$countryyear))

m1 <- function(){
  for(i in 1:N){
    propfic[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.proportional*proportional[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.proportional ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

p1 <- c("b", "b.proportional", "b.gdppercap", "b.freedomhouse")

i1 <- function(){
  list("b"=rnorm(1), "b.proportional"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}


## Criterion 1 (gallagher)
library(foreign)
dat <- read.dta("criteria.dta")
d2 <- with(dat, list(propfic=propfic, gallagher=log(gallagher), gdppercap=gdppercap, freedomhouse=freedomhouse))
d2$N <- length(unique(dat$countryyear))
m2 <- function(){
  for(i in 1:N){
    propfic[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.gallagher*gallagher[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.gallagher ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

i2 <- function(){
  list("b"=rnorm(1), "b.gallagher"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}
p2 <- c("b", "b.gallagher", "b.gdppercap", "b.freedomhouse")


## Criterion 1 (mdm)
library(foreign)
criteria <- read.dta("criteria.dta")
d3 <- with(criteria, list(propfic=propfic, mdm=log(mdm), gdppercap=gdppercap, freedomhouse=freedomhouse))
d3$N <- length(unique(criteria$countryyear))

m3 <- function(){
  for(i in 1:N){
    propfic[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.mdm*mdm[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.mdm ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
  
}

p3 <- c("b", "b.mdm", "b.gdppercap", "b.freedomhouse")


i3 <- function(){
  list("b"=rnorm(1), "b.mdm"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}

## Criterion 1 (proportional) with number of parties
criteria <- read.dta("criteria.dta")
criteria <- criteria[-which(criteria$countryyear %in% c("Thailand 2001", "Thailand 2011", "Japan 1996", "Japan 2013")),]
d4 <- with(criteria, list(propfic=propfic, proportional=proportional, gdppercap=gdppercap, freedomhouse=freedomhouse, ngov=ngov))
d4$N <- length(unique(criteria$countryyear))

m4 <- function(){
  for(i in 1:N){
    propfic[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.proportional*proportional[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]+ b.ngov*ngov[i]
  }
  b ~ dnorm(0,0.01)
  b.proportional ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  b.ngov ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

p4 <- c("b", "b.proportional", "b.gdppercap", "b.freedomhouse", "b.ngov")

i4 <- function(){
  list("b"=rnorm(1), "b.proportional"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1), "b.ngov"=rnorm(1))
}


d <- list(d1,d2,d3,d4)
p <- list(p1,p2,p3,p4)
m <- list(m1,m2,m3,m4)
ini <- list(i1,i2,i3,i4)


library(foreach)
library(doParallel)
library(R2jags)
registerDoParallel(4) 
cl <- makeCluster(4)
registerDoParallel(cl)
results <- foreach(i=1:4, .packages = c('R2jags')) %dopar% jags(data=d[[i]], model=m[[i]], parameters.to.save=p[[i]], inits=ini[[i]], n.iter=100000, n.burnin=20000, n.thin=10, n.chains=2)

results1 <- as.mcmc(results[[1]])
results2 <- as.mcmc(results[[2]])
results3 <- as.mcmc(results[[3]])
results4 <- as.mcmc(results[[4]])


# Convergence diagnostics

library(ggmcmc)
p1 <- ggs_geweke(ggs(results1)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 1")
p2 <- ggs_geweke(ggs(results2)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 2")
p3 <- ggs_geweke(ggs(results3)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 3")
p4 <- ggs_geweke(ggs(results4)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 4")
library(gridExtra)

pdf("criterion1_noties.pdf")
grid.arrange(p1,p2,p3,p4)
dev.off()

## GET OUTPUT
#Model 1 (Proportional)
coefs1 <- as.data.frame(as.matrix(results1))
b <- coefs1[,"b"]
b.proportional <- coefs1[,"b.proportional"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]

co1 <- c(mean(b), mean(b.proportional), NA, NA, mean(b.freedomhouse), mean(b.gdppercap))
sd1 <- c(sd(b), sd(b.proportional), NA, NA, sd(b.freedomhouse), sd(b.gdppercap))
p1 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.proportional)>0, mean(b.proportional>0),mean(b.proportional<0)), NA, NA, ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 2 (Gallagher)
coefs1 <- as.data.frame(as.matrix(results2))
b <- coefs1[,"b"]
b.gallagher <- coefs1[,"b.gallagher"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]


co2 <- c(mean(b), NA, mean(b.gallagher), NA, mean(b.freedomhouse), mean(b.gdppercap))
sd2 <- c(sd(b), NA, sd(b.gallagher), NA, sd(b.freedomhouse), sd(b.gdppercap))
p2 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, ifelse(mean(b.gallagher)>0, mean(b.gallagher>0),mean(b.gallagher<0)), NA, ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 3 (MDM)
coefs1 <- as.data.frame(as.matrix(results3))
b <- coefs1[,"b"]
b.mdm <- coefs1[,"b.mdm"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]

co3 <- c(mean(b), NA, NA, mean(b.mdm), mean(b.freedomhouse), mean(b.gdppercap))
sd3 <- c(sd(b), NA, NA, sd(b.mdm), sd(b.freedomhouse), sd(b.gdppercap))
p3 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, NA, ifelse(mean(b.mdm)>0, mean(b.mdm>0),mean(b.mdm<0)), ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


library(xtable)
model.effects <- data.frame(co1,sd1,p1,co2,sd2,p2,co3,sd3,p3)                                                                                                                
rownames(model.effects) <- c("Intercept", "Proportional", "Gallagher*", "MDM*", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P", "Mean", "SD", "P","Mean", "SD","P")
xtab <- xtable(model.effects, caption="Criterion 1", digits=c(2))
sink("criterion1noties.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)

#Model 4 (Proportional with ngov)
coefs1 <- as.data.frame(as.matrix(results4))
b <- coefs1[,"b"]
b.proportional <- coefs1[,"b.proportional"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]
b.ngov <- coefs1[,"b.ngov"]

co1 <- c(mean(b), mean(b.proportional), mean(b.ngov), mean(b.freedomhouse), mean(b.gdppercap))
sd1 <- c(sd(b), sd(b.proportional), sd(b.ngov), sd(b.freedomhouse), sd(b.gdppercap))
p1 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.proportional)>0, mean(b.proportional>0),mean(b.proportional<0)), ifelse(mean(b.ngov)>0, mean(b.ngov>0), mean(b.ngov<0)), ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


library(xtable)
model.effects <- data.frame(co1,sd1,p1)                                                                                                                
rownames(model.effects) <- c("Intercept", "Proportional", "Number of parties", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P")
xtab <- xtable(model.effects, caption="Criterion 1", digits=c(2))
sink("criterion1bnoties.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)


